#ifndef MODIFIERS_H
#define MODIFIERS_H
#include <stdio.h>
#include <stdlib.h>
using namespace std;

//****************************************************************************
vector<TLorentzVector> modifier_boost(long double gamma, long double theta, vector<TRootLHEFParticle*> partvector, vector<TLorentzVector> vec_partvector) {
	//****************************************************************************
	/* *********************************************  */
	/* STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ********************************************** */

	TLorentzVector jetjet;
	TLorentzVector partonePIDpositivo;

	vector<int> partons;
	partons.resize(0);

	for (int i = 1; i < 5; i++) {
		partons.push_back(i);
		partons.push_back(-i);
	}

	// bottom, top and gluon are not included


	int nmax = vec_partvector.size();

	// FIND THE INITIATING PARTONS OF THE QCD SHOWER
	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {

		int status;
		status = partvector[jcount]->Status;
		int pid;
		pid = partvector[jcount]->PID;

		if (status == -23 && isinclass(partons, pid)) {
			jetjet = jetjet + vec_partvector[jcount];

			if (pid > 0) {
				partonePIDpositivo = vec_partvector[jcount];
			}
		}
	}

	TVector3 bCoM, myboost;
	/* ******************************* */
	long double mW = 80.0;
	long double s = gamma * mW;
	/* ******************************* */

	myboost.SetX(0.0);
	// beta=sqrt(1-1/gamma^2)
	myboost.SetY(sqrt(1.0 - (mW / s) * (mW / s)) * sin(theta));
	myboost.SetZ(sqrt(1.0 - (mW / s) * (mW / s)) * cos(theta));

	bCoM = -jetjet.BoostVector(); // b=(x/t,y/t,z/t)

	//	cout << "****************" << endl;
	//
	//	print3v(bCoM);
	//	cout << "modulo di bCoM " << bCoM.Mag() << endl;
	//	cout << "modulo di myboost " << myboost.Mag() << endl;
	//	cout << "gamma di bCoM " << 1.0 / sqrt(1 - bCoM.Mag2()) << endl;
	//	cout << "gamma di myboost " << 1.0 / sqrt(1 - myboost.Mag2()) << endl;

	TLorentzVector original;
	original = jetjet;
	jetjet.Boost(bCoM);
	//	printlv(original);
	//	printlv(jetjet);

	TVector3 direction; // v2 = (1,0,0)
	TLorentzVector assez;

	assez.SetPxPyPzE(0, 0, 1, 0);
	//				cout << "partone originale";
	//				printlv(partonePIDpositivo);
	partonePIDpositivo.Boost(bCoM);
	//				cout << "partonenelcdm";
	//				printlv(partonePIDpositivo);
	direction = vectorproduct(assez, partonePIDpositivo);
	//				cout << " direction" << direction.X() << " " << direction.Y() << " " << direction.Z() << endl << "vechhio ";
	//				printlv(partonePIDpositivo);
	long double angolo;
	angolo = -partonePIDpositivo.Vect().Angle(assez.Vect());

	partonePIDpositivo.Rotate(angolo, direction);
	//				cout << "partone lungo z";
	//				printlv(partonePIDpositivo);

	//				partonePIDpositivo.Boost(bCoM);
	//				partonePIDpositivo.Rotate(angolo, direction);

	//				partonePIDpositivo.Boost(myboost);
	//				cout << "partone boostato";
	//				printlv(partonePIDpositivo);

	vector<TLorentzVector> boostedvectors;

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		TLorentzVector localLV;
		localLV = vec_partvector[jcount];
		localLV.Boost(bCoM);
		localLV.Rotate(angolo, direction);
		localLV.Boost(myboost);
		boostedvectors.push_back(localLV);
	}

	/* ****************************************************  */
	/* END OF STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ***************************************************** */

	return boostedvectors;
}

//****************************************************************************
vector<TLorentzVector> modifier_boost_relative(long double gamma, long double theta, vector<TRootLHEFParticle*> partvector, vector<TLorentzVector> vec_partvector) {
	//****************************************************************************
	/* *********************************************  */
	/* STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ********************************************** */

	TLorentzVector jetjet;
	TLorentzVector partonePIDpositivo;

	vector<int> partons;
	partons.resize(0);

	for (int i = 1; i < 5; i++) {
		partons.push_back(i);
		partons.push_back(-i);
	}

	// bottom, top and gluon are not included


	int nmax = vec_partvector.size();

	// FIND THE INITIATING PARTONS OF THE QCD SHOWER
	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {

		int status;
		status = partvector[jcount]->Status;
		int pid;
		pid = partvector[jcount]->PID;

		if (status == -23 && isinclass(partons, pid)) {
			jetjet = jetjet + vec_partvector[jcount];

			if (pid > 0) {
				partonePIDpositivo = vec_partvector[jcount];
			}
		}
	}

	TVector3 bCoM, myboost;
	/* ******************************* */
	long double mW = 80.0;
	long double s = gamma * mW;
	/* ******************************* */

	myboost.SetX(0.0);
	// beta=sqrt(1-1/gamma^2)
	myboost.SetY(sqrt(1.0 - (mW / s) * (mW / s)) * sin(theta));
	myboost.SetZ(sqrt(1.0 - (mW / s) * (mW / s)) * cos(theta));

	bCoM = -jetjet.BoostVector(); // b=(x/t,y/t,z/t)

	//	cout << "****************" << endl;
	//
	//	print3v(bCoM);
	//	cout << "modulo di bCoM " << bCoM.Mag() << endl;
	//	cout << "modulo di myboost " << myboost.Mag() << endl;
	//	cout << "gamma di bCoM " << 1.0 / sqrt(1 - bCoM.Mag2()) << endl;
	//	cout << "gamma di myboost " << 1.0 / sqrt(1 - myboost.Mag2()) << endl;

	TLorentzVector original;
	original = jetjet;
	jetjet.Boost(bCoM);
	//	printlv(original);
	//	printlv(jetjet);

	TVector3 direction; // v2 = (1,0,0)
	TLorentzVector assez;

	assez.SetPxPyPzE(0, 0, 1, 0);
	//				cout << "partone originale";
	//				printlv(partonePIDpositivo);
	partonePIDpositivo.Boost(bCoM);
	//				cout << "partonenelcdm";
	//				printlv(partonePIDpositivo);
	direction = vectorproduct(assez, partonePIDpositivo);
	//				cout << " direction" << direction.X() << " " << direction.Y() << " " << direction.Z() << endl << "vechhio ";
	//				printlv(partonePIDpositivo);
	long double angolo;
	angolo = -partonePIDpositivo.Vect().Angle(assez.Vect());

	//	partonePIDpositivo.Rotate(angolo, direction);
	//				cout << "partone lungo z";
	//				printlv(partonePIDpositivo);

	//				partonePIDpositivo.Boost(bCoM);
	//				partonePIDpositivo.Rotate(angolo, direction);

	//				partonePIDpositivo.Boost(myboost);
	//				cout << "partone boostato";
	//				printlv(partonePIDpositivo);

	vector<TLorentzVector> boostedvectors;

	myboost.Rotate(-angolo, direction);

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		TLorentzVector localLV;
		localLV = vec_partvector[jcount];
		localLV.Boost(bCoM);
		//		localLV.Rotate(angolo, direction);
		localLV.Boost(myboost);
		boostedvectors.push_back(localLV);
	}

	/* ****************************************************  */
	/* END OF STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ***************************************************** */

	return boostedvectors;
}

//****************************************************************************
vector<TLorentzVector> modifier_boost_identity(long double gamma, long double theta, vector<TRootLHEFParticle*> partvector, vector<TLorentzVector> vec_partvector) {
	//****************************************************************************
	/* *********************************************  */
	/* STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ********************************************** */

	TLorentzVector jetjet;
	TLorentzVector partonePIDpositivo;

	vector<int> partons;
	partons.resize(0);

	for (int i = 1; i < 5; i++) {
		partons.push_back(i);
		partons.push_back(-i);
	}

	// bottom, top and gluon are not included


	int nmax = vec_partvector.size();

	// FIND THE INITIATING PARTONS OF THE QCD SHOWER
	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {

		int status;
		status = partvector[jcount]->Status;
		int pid;
		pid = partvector[jcount]->PID;

		if (status == -23 && isinclass(partons, pid)) {
			jetjet = jetjet + vec_partvector[jcount];

			if (pid > 0) {
				partonePIDpositivo = vec_partvector[jcount];
			}
		}
	}

	TVector3 bCoM, myboost;
	/* ******************************* */
	long double mW = 80.0;
	long double s = gamma * mW;
	/* ******************************* */

	//	myboost.SetX(0.0);
	//	// beta=sqrt(1-1/gamma^2)
	//	myboost.SetY(sqrt(1.0 - (mW / s) * (mW / s)) * sin(theta));
	//	myboost.SetZ(sqrt(1.0 - (mW / s) * (mW / s)) * cos(theta));

	bCoM = -jetjet.BoostVector(); // b=(x/t,y/t,z/t)

	//	cout << "****************" << endl;
	//
	//	print3v(bCoM);
	//	cout << "modulo di bCoM " << bCoM.Mag() << endl;
	//	cout << "modulo di myboost " << myboost.Mag() << endl;
	//	cout << "gamma di bCoM " << 1.0 / sqrt(1 - bCoM.Mag2()) << endl;
	//	cout << "gamma di myboost " << 1.0 / sqrt(1 - myboost.Mag2()) << endl;

	TLorentzVector original;
	original = jetjet;
	jetjet.Boost(bCoM);
	//	printlv(original);
	//	printlv(jetjet);

	TVector3 direction; // v2 = (1,0,0)
	TLorentzVector assez;

	assez.SetPxPyPzE(0, 0, 1, 0);
	//				cout << "partone originale";
	//				printlv(partonePIDpositivo);

	long double angoloLAB;
	angoloLAB = partonePIDpositivo.Vect().Angle(original.Vect());

	long double pTrelpidpositivo;
	pTrelpidpositivo = partonePIDpositivo.Vect().Mag() * sin(angoloLAB);

	long double senThetaCoM;
	senThetaCoM = pTrelpidpositivo / (original.Mag() / 2.0);

	/* ****************** */
	/* BOOSTO IL PARTONE COL PID POSITIVO */
	partonePIDpositivo.Boost(bCoM);
	/* ****************** */

	// angolo di ROOT nel CoM
	long double angoloROOT;
	angoloROOT = partonePIDpositivo.Vect().Angle(original.Vect());

	/* ********************* */
//	theta = angoloROOT;
	/* ********************* */
//		cout << "************" << endl;
//	cout << angoloROOT << endl;
//	cout << theta << endl;

	myboost.SetX(0.0);
	// beta=sqrt(1-1/gamma^2)
	myboost.SetY(sqrt(1.0 - (mW / s) * (mW / s)) * sin(theta));
	myboost.SetZ(sqrt(1.0 - (mW / s) * (mW / s)) * cos(theta));


	//	cout << " angolo di ROOT nel CoM " << angoloROOT << endl;
	//	cout << " sin angolo nel CoM calcolato dal LAB " << senThetaCoM << endl;
	//	cout << " massa"<<endl;

	//				cout << "partonenelcdm";
	//				printlv(partonePIDpositivo);
	direction = vectorproduct(assez, partonePIDpositivo);
	//				cout << " direction" << direction.X() << " " << direction.Y() << " " << direction.Z() << endl << "vechhio ";
	//				printlv(partonePIDpositivo);
	long double angolo;
	angolo = -partonePIDpositivo.Vect().Angle(assez.Vect());

	myboost.Rotate(-angolo, direction);


//	partonePIDpositivo.Rotate(angolo, direction);
	//				cout << "partone lungo z";
	//				printlv(partonePIDpositivo);

	//				partonePIDpositivo.Boost(bCoM);
	//				partonePIDpositivo.Rotate(angolo, direction);

	//				partonePIDpositivo.Boost(myboost);
	//				cout << "partone boostato";
	//				printlv(partonePIDpositivo);

	vector<TLorentzVector> boostedvectors;

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		TLorentzVector localLV;
		localLV = vec_partvector[jcount];
		localLV.Boost(bCoM);
//		localLV.Rotate(angolo, direction);
		localLV.Boost(myboost);
		boostedvectors.push_back(localLV);
//		cout << "vettori " << jcount << "esimi" << endl;
//		printlv(vec_partvector[jcount]);
//		printlv(boostedvectors[jcount]);
	}

	/* ****************************************************  */
	/* END OF STUFF TO BOOST THE PARTICLES BEFORE CLUSTERING */
	/* ***************************************************** */

	return boostedvectors;
}

#endif
